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Abstract 

(— ( Impulse methods are generalized to a family of integrators for Langevin systems with 

r-{ quadratic stiff potentials and arbitrary soft potentials. Uniform error bounds (independent 

from stiff parameters) are obtained on integrated positions allowing for coarse integration 
steps. The resulting integrators are explicit and structure preserving (quasi-symplectic for 
Langevin systems). 

1 Introduction 

Results: This paper generalizes the impulse methods for stiff Hamiltonian systems [TSJ [32] to 
stiff stochastic Langevin systems. In the stochastic setting these integrators are quasi-symplectic 
as defined in [29]. Unform error bounds are obtained for both stochastic and deterministic 
settings. 

More precisely, this paper is concerned with the numerical integration the following stiff 
SDEs: 

f Mdq = pdt 

^ \ dp = -VV(q)dt-e- 1 Kqdt~cpdt + crdW 

which describes a stochastic mechanical system with a potential being sum of slow V(q) and fast 
\e~ x q T Kq, and the momentum being perturbed by noise and attenuated by friction. 

When noise and friction are both present (i.e. c ^ 0, a ^ 0), a lst-order member of the 
proposed Stochastic Impulse Methods (SIMs) family will integrate position q with a global error 
uniformly bounded by CH 1 / 2 , where H is the integration timestep and C is a constant indepen- 

y—{ dent from e _1 (provided that the solution remains bounded). The integrator is also shown to be 

quasi-symplectic. When noise and friction are absent (i.e. c = 0, a = 0), the deterministic (and 

• i-H symplectic) version of the lst-order SIM gives a uniform lst-order global error on q, if again the 

solution is bounded. The error bound on momentum p, however, is not uniform here. Recall that 
increased accuracy and stability was one of motivations supporting the development of mollified 
impulse methods [TU [32] . 

Dynamical systems with multiple time scales pose a major problem in simulations because the 
small time steps required for stable integration of the fast motions lead to large numbers of time 
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steps required for the observation of slow degrees of freedom |16| . As seen from the error bounds, 
in the case of quadratic fast potential, SIMs provide a possibility of accurate integration with 
a choice of timestep not restricted by the stiffness e _1 , as long as position is the quantity of 
interest. In these large timestep can be adopted. 

Also, SIMs are symplectic [TB] and in fact variational in the case of no noise no fric- 
tion, and are quasi-symplectic [35J in the case of full Langevin. As a result of the preservation of 
structure, properties such as near preservation of energy or of associated Boltzmann-Gibbs invari- 
ant measure, as well as conservation of momentum maps could be obtained, which significantly 
benefit long time numerical integrations. 

Related work: Many elegant methods have been proposed in the area of stiff Hamiltonian/Langevin 
integration with different focus and perspective. 

Impulse methods, as well as other members of the exponential integrator family [T3] , including 
Mollified Impulse Methods [12] [32] and Gautschi-type integrators [17] are prevailing symplectic 
integrators for stiff Hamiltonian systems. They are however not directly extendable to stiff 
Langevin systems if integration with a large step is desired. The general GLA [3J approach 
(see also [31] for an extension of impulse methods non-stiff Langevin systems) of constructing 
Langevin integrator from a symplectic scheme by composing an Ornstein-Uhlenbeck flow with 
the symplectic integrator will not yield a uniform error bound in the case of stiff Langevins. It 
is worth mentioning that impulse methods are not limited to quadratic stiff potentials (provided 
that the flow of the stiff part of the Hamiltonian is given). 

The implicit method approach for integrating stiff equations include the LIN algorithm [ID] 
for stiff Langevin systems. However, it has been observed that Implicit methods in general fail to 
capture the effective dynamics of the slow time scale because they cannot correctly capture non- 
Dirac invariant distributions |24j . Moreover implicit methods are generally slower than explicit 
methods, provided they use comparable timesteps. 

Implicit and explicit approaches were combined in a variational integration framework by 
defining the discrete Lagrangian via trapezoidal approximation of the soft potential and midpoint 
approximation of the stiff potential. The resulting IMEX for stiff Hamiltonian systems [35] is 
explicit in the case of quadratic fast potential. Similar as the case of impulse methods, there is 
no easy way to extend IMEX to stiff Langevin systems. 

The Hamilton-Jacobi derived homogenization method for multiscale Hamiltonian systems |22j 
enables the usage of a large timestep for deterministic systems but can not directly be extended 
to stiff Langevin systems using the GLA approach [3J . 

Multiscale methods that integrate the slow dynamics by averaging the effective contribution 
of the fast dynamics have been applied to stiff Langevin systems. These include Heterogeneous 
Multiscale Methods (HMM) [TOl HJ [1] , equation free methods [HI H3j [20] , and FLow AVer- 
aging integratORS (FLAVORS) [37]. We observe that these methods use mesoscopic timesteps, 
which are (usually) one or two orders of magnitude smaller than the large steps employed in 
SIMs. These methods also assume a separation of timescales, and therefore will not work for 
generic stiff Langevin systems that are not necessarily multiscale. In addition, based on averag- 
ing instantaneous drifts, both Heterogeneous Multiscale Methods and equation free methods (in 
their original form) require an identification of slow variables in general nonlinear cases (with 
exceptions such as in (DJ). Reversible and symmetric methods in these frameworks have been 
proposed [53J [U [33J . FLAVORS are based on averaging instantaneous flows and do not require 
explicit identification of slow variables, and are symplectic (quasi-symplectic). 
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2 Stochastic Impulse Methods 



Consider numerical integration of the following multiscale Langevin SDEs 



Mdq = pdt 

dp = -W(q)dt- e^Kqdt- cpdt + adW 



(1) 



where < e <C 1, q E M. d , p E R d , K is positive definite d x d matrix, c and a are positive 
semi-definite d x d matrices, respectively indicating viscous damping coefficients and amplitudes 
of noises. We restrict ourselves to Euclidean phase spaces, although the method is readily 
generalizable to manifolds. In addition, we require that matrices K and c commute; a special 
case satisfying this requirement is c being a scalar. 

In the case of no noise no friction (c = and a = 0), the system degenerates to a deterministic 
mechanical system with Hamiltonian H(q,p) — ^p T M~ 1 p + V(q) + eT x \q T Kq. 

Also, the method as well as the uniform convergence theorem works for a more general open 
system: 

Mdq = pdt 

dp = F(q)dt - e- 1 Kqdt - cpdt + adW 

but we stick to ^ for simplicity in descriptions. 

Denote by <^(r) : (qf \t),pf '(*)) ^ (q* (t + r) , pf (t + r)) and0 s (r) : (q s (t) , p s (t)) ^ (<f(i + r),p s (i + r)) 
respectively the r-flow maps of the autonomous SDE systems 



(2) 



Mdq f 
dpf 



= pfdt 

= -e- 1 Kqi'dt-cp f dt + adW 



and 



Mdq s 
dp s 







-VV{q s )dt 



(3) 



(4) 



Since the first system is a linear SDE and the second is a free drift, flows of both can be obtained 
exactly. 

Then Stochastic Impulse Methods(SIMs) are defined via compositions of <f>* and <fi 3 . Here 
are several examples of SIMs with a timestep H: 



Integrator 1. 1st order SIM in the c 

(f> s (H) o (f>f(H): 



0, a = case, is given by the one step update of 



where 



Pk' 

Qk+l 

Pk+l 

A U (H) A 12 (H) 
A 21 (H) A 22 (H) 



A n (H)q k + Ai 2 (H)p k 
A 2 i{H)q k + A 22 (H) Pk 
qk' 

p k ,-HVV{q k .) 

\ 

= exp _ e -i KH 



M- X B 




I 9o = q(0) 
\po=p(0) 



Remark 2.1. The other 1st order SIM, as the above's dual, can be obtained via the one step 
update & {H) o <f> s (H). Both these 1st order composition schemes are well known as the Lie- 
Trotter splitting [38]. 
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Integrator 2. 1st order SIM in the full Langevin case, given by the same one step update 
(f> s (H) o (f> s (H): 



qk' 
Pk> 

Qk+l 

Pk+i 
Rq k (H) 
Rpk(H) 



where 



Bu{H)q k 
B 2l {H)q k 
qk' 

py - HVV(q k 





B 12 {H)p k 
B 22 {H)p k 



Rq k (H) 
Rpk(H) 



E2,(£f)) 



B n (H) 
B 2 i(H) 



2 

21 

B 12 (H) 
B 22 (H) 



Y 2 

^22 



(H)) 



),i.i.d. normal distributed 



exp 







M~ X H 

-cH 



Efi(-ff) 
£?,(ff) 



I 9o = 9(0) 
\po=p(0) ' 

J" (B 12 (H-s)oa T BUH 



s)) ds 



Co {B 12 (H-s)aa T BUH~s)) ds 
(H) - X,= (B 22 (H - s)aa T Bf 2 (H s)) ds 
%h(H) = J s t (B 22 (H - s)aa T BUH ~ a)) ds 



J 12V 

Y 2 



Remark 2.2. 



Rq k (H) 

RPk(H) 



indicates the value of J s=0 B(H — s 



normal random variable with zero mean and covariance of 





adW. 
YUH) Y\ 2 (H 

•2 
22 



and hence is a vectorial 



Y 2 21 (H) Y 2 22 (H) 



Integrator 3. 2nd order SIM in the full Langevin case, given by the one step update 
(j> s (H/2) o (f>f(H) o (f> s (H/2): 



qk' = qk 

Pk> = Pfc-fW(g fe ) 

q kl , = B n (H)q k , + B 12 (H)p k , 

Pk>> = B 2l {H)q k , + B 22 {H)p k , 

qk+i = qk" 

Pk+i = Pk" - f-W(g fc //) 



Rq k {H) 
Rpk(H) 



Remark 2.3. This uses the 2nd order composition scheme known as the Strang or Marchuk 
splitting [23]. When no noise or friction, i.e. c = and X = 0, the resulting integrator 
degenerates to the prevailing Verlet-I/r-RESPA impulse method [T51 139) . 

Remark 2.4. Higher order SIMs can be obtained systematically since generic way for construct- 
ing higher order splitting/composition schemes exists [IB]. For instance a 4th order SIM is given 
by (f> s (cH/2) o (f)f(cH) o </> s ((l - c)H/2) o </>f((l - 2c) H) o S ((1 - c)H/2) o ^ (cH) o (/) s (cH/2) 
where c = 2 _l 1/3 [3D] . 



3 Properties 
3.1 Symplecticity 

In the case of c — and a = 0, since <j) s and $ are the exact flows of Hamiltonian systems, they 
are symplectic. Therefore SIMs, as compositions of the two, are symplectic. 
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In fact, SIMs here are not only symplectic but variational, in the sense that their equations of 
motion are obtained as critical point of a globally defined action, which is the integral of a discrete 
Lagrangian. Since SIMs are based on splitting schemes, and the original system is split to two 
Hamiltonian systems, backward error analysis can be done via Poisson brackets [16] . resulting in 
a global non-dcgcncratc Hamiltonian that is exactly preserved. Then Legendre transformation 
gives the discrete Lagrangian and hence the variational structure. 

When noise and friction are present, SIMs are quasi-symplectic for RL1 and RL2 in [33] can 
be easily checked to be true, i.e. they degenerate to symplectic methods if friction is set equal 
to zero and the Jacobian of the flow map is independent of (q,p). 

If in addition c is isotropic, then SIMs are conformally symplectic, i.e. they preserve the 
precise symplectic area change associated to the flow of inertial Langevin processes [25] . 



3.2 Uniform Convergence 

In the case of c = and a = 0, convergence of SIMs is guaranteed by the general construction 
of splitting schemes. In the full Langevin setting, analogous convergence results for the same 
splitting schemes can be easily obtained using generators of SDEs. By this approach, however, 
the error bound will contain the scaling factor e~ 1 and therefore restrain the timestep from being 
large. We instead seek for uniform convergence results, i.e. error bounds that don't depend on 
Li. It turns out such a uniform bound holds only for the position q but not the momentum p. 

Condition 3.1. We will prove a uniform bound on the scaled energy norm of the global error 
of Integrator^ if the following conditions hold: 

1. Matrices c and K commute. A special case could be c being a scalar. 

2. lim^o -y/e 1 1 c| 1 2 < C for some constant C independent of e, i.e. c < 0(e -1 ^ 2 ). 

3. a is independent of e _1 , in the sense that lim e _^.o e p 1 1 cr 1 1 2 = for any p > 0. 

4- In the integration domain of interest VV(-) is bounded and Lipschitz continuous with coef- 
ficient L, i.e. || W(a) - W(6)|| 2 < L\\a - b\\ 2 . 

5. Denote by x(T) — (q(T),p(T)) the exact solution to and xt — (<?t,Pt) the discrete 

I^tIH — C f or some 



numerical trajectory given by Integrator^ then ¥,\\x(T)\\2 < C and E 
constant C independent of e _1 but dependent on initial condition E| 



2, amplitude of 
noise a and friction c. 

Note that this condition usually holds due to preservation of Boltzmann-Gibbs invariant 
measure, whose parameter of temperature doesn't depend on e _1 since noise is weak, and 
whose energy function is usually dominated by the positive definite fast potential ( implying 
closed energy level sets). 



Theorem 3.1. If Condition 3.1 holds, the 1 order SIM (Integrator^ for multiscale Langevin 
system ^ (c ^ 0, a ^ 0) has in mean square sense a uniform global error of 0(i? 1//2 ) in q and 
a non-uniform global error of €~ 1 / 2 0{H 1 / 2 ) in p, given a fixed total simulation time T — NH: 

(EH^-gTll 2 ) 1 / 2 < CH 1 ' 2 (5) 
(EUp^-ptH 2 ) 1 / 2 < e-^WVKhCH 1 / 2 (6) 



where q(T),p(T) is the exact solution and qT,PT is the numerical solution; C is a positive 
constant independent of t~ x but dependent on simulation time T, scaleless elasticity matrix K, 
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scaled damping coefficient ^fic (0(1)), amplitude of noise a, slow potential energy V(-), and 



initial condition Ell 



Proof. We refer to the appendix for the proof. 



□ 



Remark 3.1. By looking at the proof, one can be assured that all convergence results of 

,2 

SIMs apply to situations where the deterministic system is in a more general form of M-^q — 
—e~ 1 Kq + F(q), where F(q) doesn't have to be — W(<?). 

In the special case of Hamiltonian system, the same integrator gains 1/2 more order of 
accuracies. 

Condition 3.2. We will prove a uniform bound on the scaled energy norm of the global error 
of Integrator^ if the following conditions hold: 

1. In the integration domain of interest VV(-) is bounded and Lipschitz continuous with coef- 
ficient L, i.e. || W(a) - W(6)|| 2 < L\\a - b\\ 2 . 

2. Denote by x(T) = (q(T),p(T)) the exact solution to ^ with c = and a — 0, and 
X T = (qt,Pt) the discrete numerical trajectory given by Integrator^ then \\x(T)\\2 < C 

£t||! < C for some constant C independent of e _1 but dependent on initial condition 



and 
Pa 



Note that this condition usually holds due to preservation of energy, which is usually dom- 
inated by the positive definite fast potential (implying closed energy level sets). 

Theorem 3.2. If Condition 3.2 holds, the 1 st order SIM (Integrator^ for multiscale Hamilto- 
nian system ( ^ with c — 0, a = 0) has a uniform global error of O(H) in q and a non-uniform 
global error of e" 1 / 2 0(H) in p, given a fixed total simulation time T — NH: 



UT) 
\P(T) 



■ qr\\2 
Prh 



< 
< 



CH 



Vk\uch 



(7) 
(8) 



where q(T),p(T) is the exact solution and qT-,PT is the numerical solution; C is a positive 
constant independent of e _1 but dependent on simulation time T, scaleless elasticity matrix K, 

slow potential energy V(-) and initial condition \\ ^° 
Proof. It follows by simplifying the proof of Theorem |3.1[ □ 



3.3 Stability 



As one sees from Condition 3.1 and 3.2 (as another nonlinear demonstration of Lax equivalence 
theorem [21) ). stability is necessary for global convergence. Instability could either come from 
the problem itself (not all SDEs have bounded solutions in the mean square sense), or from 
imperfection in numerical integration schemes. Here consider the latter possibility only. It is 
shown that impulse methods are not unconditionally stable [12] . and its improvement, mollified 
impulse methods, are still susceptible to instability intervals (although narrower) in a linear 
example [5]. Nevertheless, instability intervals of impulse method are already narrow regions; for 
instance, the first instability interval in the stiff example considered by [5] is 0.544 < H < 0.553. 
It is intuitive that instability intervals for the stochastic case with damping or higher order 
schemes will not be wider. Therefore one could still choose a large timestep H in SIMs without 
hitting the instability, by at most a few integration tryouts with slightly varied i?'s. 



G 



4 Numerical Examples 

4.1 2-spring systems with noise and friction 




Figure 1 : 2-spring systems 

Consider a "Wall - linear stiff Spring - Mass - nonlinear soft Spring - Mass" system with 
both masses under isotropic noise and friction (Figure 0. The governing equations write as: 

dx — p x dt 

dy — pydt 

dp x = —(uj 2 x + (x — y) 3 )dt — cp x dt + adW} 

k dp y — — (y — x) 3 dt — cpydt + adWj 2 

Note (1) this is a Langevin system with H(x, y,p x ,p y ) — \p x + \p 2 y + \ijPx 1 + \{y — x) 4 (2) 
y is a slow variable but x is not purely fast (there is a slow component in it). 




Time Time 



(a) FuJJ period case: s'm(uiH) = (b) Quarter period case: cos(ujH) = 

Figure 2: Empirical moments obtained by lst-order SIM with large step H and lst-order GLA [3] with small 
step h. Parameters are w = 100, c = 0.1, P = |§ = 10, x(0) = 0.8/oj, y(0) = 1.1 + x(0), p x (0) = 0, p y {0) = 0; 
h = 0.1/lj and H is chosen to be not scaling with ui yet corresponding to a resonant frequency; empirical moments 
are obtained by averaging 5000 simulations. 

lst-order SIM (Integrator [2]) is compared in Figure[2]to the benchmark of Geometric Langevin 
Integrator (GLA) [3 which is Boltzmann-Gibbs preserving and convergent. Agreements on 
empirical moments of integrated trajectories serve as evidences of structure preservation and 
convergence. The large timesteps used by SIM are chosen to be the resonance frequencies and 
they do produce stable accurate results. 0(w)-fold acceleration is gained by SIM. 
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it la • • • q 2 ».i q 2m 

stiff soft I 

Figure 3: Fermi-Pasta-Ulam problem 11 - ID chain of alternatively connected harmonic stiff and non-harmonic 
soft springs 

4.2 Fermi-Pasta-Ulam problem 

Consider the deterministic Fermi-Pasta-Ulam (FPU) problem [11] illustrated in Figure [3] and 
associated with the Hamiltonian 

h ^p) - 2^ ( ^-i+^) + ^I](fe-fe-i) 2 + ^(g2 l +i-fe) 4 (9) 

1 = 1 1 = 1 2=0 

Conventionally the following transformation is used 

aft = (fe + 92i-i)/v / 2 

Xm+i = (fe-fe-i)/v| i = l m (10) 

Hi = (P2i +p 2 i-l)/v2 
2/m+J = (P2» -P2i-l)/v / 2 

so that the fast potential is diagonalized: 

'H(x,y) =\Y*™xyt + V f {x)+V s {x) 



2 

m-\-i 



v f (x) =\Y™i 

V s (x) = \{{xi - X m+1 ) 4 + YX=l 1 ( X i+l ~ x m+i+l - Xi ~ X m+i ) A + (x m + X 2m ) 4 ) 

The FPU problem is a well known benchmark problem |271 116) for multiscale integrators 
because it exhibits different behaviors over widely separated timescales. The stiff springs (nearly) 
behave like harmonic oscillator with period ~ 0(w^ 1 ). Then the centres of masses linked by 
stiff springs (i.e. the middle points of stiff springs) change over a timescale 0(1). The third 
timescale 0(lj) is associated with the rate of energy exchange among stiff springs. On the other 
hand, in addition to conservation of energy, the total energy of stiff springs behave almost like 
a constant. Comprehensive surveys on FPU problem, including discussions on timescales and 
numerical recipes, can be found in [161 [5]. 

We present in Figure [4] lst-order SIM simulation (Integrator [lj together with variational 
Euler (a.k.a. symplectic Euler) simulation of FPU over a time span of 0(uj). Good results are 



obtained by SIM beyond the timescale of O(l) (as guaranteed by Theorem 3.1 ) but actually over 
0(ui), and 200-fold (to = 200) acceleration is gained at the same time. 

Notice that Mollified Impulse Methods with ShortAverage, LongAverage or LinearAverage 
filters [T2] didn't accurately capture the rates of energy exchanging among stiff springs over 
T = 0(uj) (results not shown). 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.9 1 
Middle point position of the first stiff spring 




5 10 15 20 25 30 35 40 45 50 
Energies of stiff springs 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.3 0.9 1 
Middle point position of the first stiff spring 




5 10 15 20 25 30 35 40 45 50 
Energies of stiff springs 




100 200 300 400 500 BOO 700 800 900 1000 100 200 300 400 500 BOO 700 800 900 1000 



oia! energy of stiff springs 



otel energy of stiff springs 



100 200 300 400 500 BOO 700 800 900 1000 [I 100 200 300 400 500 BOO 700 800 900 1000 



(a) lst-order SIM, large step H = 0.1 



(b) Variational Euler, small step h = 0.1/cj 
0.0005 



Figure 4: Simulations of FPU over T = 5w. Parameters are w = 200, m = 3, x(0) = [1,0,0, l/w,0,Q], 
y(0) = [0,0,0,0,0,0]. Different subplots use different time axes to accentuate different timescales: Subplotl 
shows scaled expansions of three stiff springs x m +i, which are fast variables; Subplot2 shows scaled middle point 
position of the first stiff spring xi, which is one of the slow variables; Subplot3 shows the energy transferring 
pattern among stiff springs, which is even slower; Subplot4 shows the near-constant total energy of three stiff 
springs. The fast variables of stiff spring expansions are in fact oscillating much faster than shown in Subplots 1, 
for Subplots 1 are plotted by interpolating mesh points with a coarse mesh size of H. 



6 Appendix 

6.1 Proof of Theorem 13.11 

Throughout this subsection Condition |3.1| is assumed. For a concise writing we also abuse the 
notation 0(x n ), which indicates some entity whose norm < Cx n , where C is a constant that 
doesn't change with e, i.e. doesn't depend on e -1 . 

Definition 6.1. Scaled energy norm: 



q 




q 


P. 


IU = II 





h = \/l T< l + ep T K~ 1 p 



This is well defined because K is positive definite. 

Since e is very small, the following inequalities for converting between scaled energy norm 
and two-norm can be easily obtained: 



Proposition 6.1. Let x = 



be any vector, then 



e 1 / 2 ||v / ^|| 2 - 1 IN| 2 = IN2 1 INl2<INU< ||z||a 



< HIT 1 ! 



= e 1/2 \\VK 



A" 'ihlklU 



Also, vector-norm-induced matrix norms satisfy 

\\Mx\\ E 



Mn Mia 
M 2 i M 22 



E = SUp 



Me 



Mh Miafi 

Q^Mai n- 1 M 22 n 



(11) 
(12) 

(13) 
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Lemma 6.1. Let B(s) = 



B u (s) B 12 (s) 
B 21 {s) B 22 {s) 



exp(s 







defined in Integrator^ with H — s and arbitrary k, then 



), and Rq(s) be the Rqk(H) 



\\Bn{s)h 


< 


1 


(14) 


ll-B 2 2(s)||2 


< 


1 


(15) 


l|Bia(*)||a 


< 


N 


(16) 


4B 2 i{s)h 


< 


C K \s\ 


(17) 


e^WBnis) - I\\ 2 


< 


C k \s\ 


(18) 


e 1/2 \\B 22 (s) - I\\ 2 


< 


C c \s\ 


(19) 


E\\Rq(s)\\ 2 2 


< 


\\W\\l\s? 


(20) 


e 1/2 \\B(s)-I\\ 2 


< 


C Kc \s\ 


(21) 


e^\\B{s)-I\\ E 


< 


CkM 


(22) 



where Ck, C c and Ckc o,re some positive real constants (may indicate different values in different 
inequalities), respectively dependent on K, \ftc, K and yfec but independent of e _1 . 

Proof. Since c and K commute, they can be diagonalized simultaneously |18j . By the theory of 
linear ordinary differential equations [31 , one can hence diagonalize B Ul B± 2 , B 2 i, B 22 simul- 
taneously. Since each diagonal element can be investigated individually, assume without loss of 
generality that = = e~ x l 2 \fK and c are both scalars, and use the notation of scalar w 
and scalar c thereafter. 

Denote the damping ratio by £ = tj. The solution to damped harmonic oscillator can be 
analytically obtained, and hence components of the flow operator Bu,B 12 ,B 2 \,B 22 as well. 
When C < 1 i- c - undcrdamping, which is usually the case since w is large 



Bu(s) 

B 12 (s) 

B 21 (s) 
B 22 (s) 



^ s (cos(u J ^/l-C 2 s) + ^£ 
V 1 _ 



=«n(Vl - C 2 *)) 



e u "' s sin 



(uj^T^Cs) 



C 2 s) 



e-^ s (cos{cj^l-C 2 s) 
When C = 1 i- e - critical damping, 



zsin(\/l — ( 2 s)) 



Bn(s) 
B 12 (s) 
B 21 (s) 
B 22 (s) 



= e- us (l+Lus) 



! (1 



(23) 

(24) 

(25) 
(26) 



(27) 
(28) 
(29) 
(30) 
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When £ > 1 i.e. over damping, 

A(s) = e ^(-C-V<^) (31) 
B(s) 4 e ««(-C+v^I) (32) 



»„(.) = « B ^» 2 + J1?^ + B » (33, 

B 21 (s) = \ ; 35 
2^(2^ 

= (36) 

(37) 

By routine investigations on local extremes using calculus, it can be shown in all three cases 
that 



\\Bu(b)\ 


2 


< 


1 


\\B 2 2(S)\ 


2 


< 


1 


\\Bu(s)\ 


2 


< 


S 


\\B2i(s)\ 


2 


< 


u 2 s 


Bn(8)-I\ 


2 


< 


US 


B 2 2(s)-I\ 


2 


< 


j US 

[ 20 



C< i 



(38) 
(39) 
(40) 
(41) 
(42) 

(43) 



When C > 1, since c = O^ 1 / 2 ) (Condition ^ , 2(us = 0{e~ l / 2 )s. Therefore e 1/2 ||£ 22 (s) 
I\\2 < C c \s\ always holds. 
Also, 



E\\Rq(s) 



Eli 



B l2 {t)adW t \\i 
= J'\\<rB 12 mlto<\Ml\*\ 



(44) 



For a proof on norm bounds of the entire matrice we use only bounds of dimensionless block 
elements: 



\B - 1\ 



< 



Bu — I B12 
B21 B22 — I 

n 



n 

1/2 



n-^Bn-I) Q,~ l B 12 



O(s) 
e-V20(« 



n- x B 21 n-\B22-i) 

e l ' 2 0{s) 



0(a) 



It's easy to prove that for any scalar a 





aMx2 




'M n 


M 12 ~ 




M21 


M22 _ 


h = \\ 


aM 2 \ 


M 22 _ 


II2 



(45) 
(46) 
(47) 

(48) 
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Therefore 



Similarly, 



\B-I\\ 2 = e-^\ 



\B - I\\ E = e- 1/2 | 
< e- 1/2 | 



0(a) 0(a) 
0(a) 0(a) 



; 2 = e 



(49) 



= e -V2 



n- 1 B 21 n- 1 B 22 n-i 

0(s) e 1 /2 e -V2o( s y 

e -i/2 e i/2 0(s) e i/2 0(s)e -i/2 
0(a) o( s y 



0(s) 0(s) 



(50) 

□ 



Remark 6.1. In the special case of c = 0, bounds of block elements can be easily obtained since 

|cos(ws)| < 1 



eK 1 | — tusin(ujs)\ 



-uain(ujs) 



< 8 



e 1 / 2 \/K \C08 (ojs) - 1| = | - 2sin 2 (ujs/2)/uj\ < \ - 2ain(oja/2)/oj\ < \s\ 

Lemma 6.2. The solution to SDE dX = AXdt + f(X)dt + HdWt can be written in the following 
integral form: 

X(t) = e At X(0) + f e A ^- s \f(X{ S ))ds+ f e A(t - s) ZdW s 
Jo Jo 



(51) 



Proof. Let Y(t) = e - At X(t), then by Ito's formula and dX = AXdt + f(X)dt + EdW t 



dY = e- At f(X(t))dt + e- At Y,dWi 



At s 



This in the integral form is 

Y(t)=Y(0)+ f e- As f(X(s))ds+ f e~ As T,dW s 
Jo Jo 



Hence 



X(t) = e At X(0) + f e A{t - s \f(X(s))ds+ f e A(t ~ s ^dW s 
Jo Jo 



(52) 

(53) 

(54) 
□ 



Lemma 6.3. Consider two continuous stochastic dynamical systems, the original dynamics and 
the bridge dynamics: 



(55) 



dq 


= pdt 


dp 


= 


5(0) 


= qo 


I P(0) 


= Po 
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dq = pdt 

dp = -e- 1 Kqdt-W{q a )dt-cpdt + adW t 

?(0) = go 

P(0) = Po 



(56) 



TTiera (E|| 



p(iO 



Hi,) 1 / 2 < C\H\ 3 / 2 , where C is a positive constant independent of e 1 



but dependent on the scaleless elasticity matrix K , scaled damping coefficient y/ec, amplitude of 
noise a, and slow potential V(-). 



Proof. Rewrite the original dynamics ( 55 ) as 



dq ~ pdt 

dp = e- 1 Kqdt-VV{q )dt+(yV(q )-VV{q))dt-cpdt + adW t 

<Z(0) = qo 

{ P(0) = Po 



Let x(t) 



, x{t) 



'«(*) 
p(t) 
o 

W(? ) ~ VV{q) 



'«(*)' 
P(*). 

and £ = 



, E(i) = exp(t 



/ 



Then by Lemma 



and bridge dynamics can be respectively written as: 



6.2 



), 6 = 





-W(<7o) 



(57) 



, fffop) = 



solutions to the original dynamics 



x(t) = B{t)x(0)+ [ B(t-s)bds+ [ B(t- s)XdW s + [ B{t - s)g(x(s))ds 
Jo Jo Jo 

x(t) = B(t)x(0)+ [ B(t-s)bds+ [ B{t-s)Y,dW s (58) 
Jo Jo 

Notice for any vector y and positive t that ||-B(t)y||E < \\v\\e, because energy is decaying in 
the system q + cq + e~ 1 Kq — 0. Together with Cauchy-Schwarz we have 



E\\x(t) -x(t)\ 



El 



B(t- s)g(x(s))ds\ 



< t I E\\B(t - s)g(x(s))\\%ds 



< t / E\\g(x(s))\\ E d S 
Jo 



(59) 



By Condition 3.1 assume VV(-) is Lipschitz continuous with coefficient L, then almost surely 

\\g(x(s))\\ E = 



o 



VV(q ) - VV(q(s)) 



< ^\\VK h\\VV(qo)-VV(q(s))\\r 



< L^~e\\VK \\ 2 \\q(s) - q \\ 2 



Similarly, since 



x(t) = B(t)x(0) + / B(t-s) 
Jo 





■vv(?(»)) 



ds+ [ B(t-s)Y,dW s 
Jo 



(60) 



(61) 
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we have 



\x(s) - B(s)x 



f B(s - t)VdWth < I \\VV{q(t))\\ 2 dt 
Jo Jo 



(62) 



By Condition 3.1 VV(-) is bounded, and hence the above is O(s) 



We now can bound (60) and therefore (59) with the aid of ( |62[ ) and Lemma 6.1 

E|k(«)-9bHa 
< E||a;(s) -xolla 



< E (\\x(s) - B(s)x - B(t 



s)XdW s \\ 2 + \\B(s)x Q - x || 2 + || / B{t-s)ZdW s \ 



< 3E ( ||x(s) - B(s)x - / B(t -s)SdW s || + ||B(s)x -x ||2 + || / B(t-s)SdW s || 



3^O(s")+e- i O(^)E||x ||^+ J a z (B 12 (t-sy +B 22 (t-sY)dt 
0(s 2 ) + e- 1 O(s 2 )E||x || + 0(s 3 ) + O(s) 



(63) 



By Condition 3.1 E||xo||2 = 0(1). Therefore, the above expression is e 1 0(s 2 ) + 0(s). 
This gives E|j(?(x(s))|j!j — O(s) independent of e 1 , and eventually E||x(/i) — x(K)\\\ — 
0{h 3 ). ' □ 

Lemma 6.4. Consider the discrete stochastic dynamical system given by Ist-order SIM (Inte- 
grator^: 



Ph 



B n (H)q + B 12 (H) Po + Rq{H) 

B 21 (H)q + B 22 (H) Po + Rp(H) - HVV(B n (H)q + B 12 (H) Po + Rq(H)) 



(64) 



Then a comparison with bridge dynamics ( |56[ ) gives E\\qn — q~{H)\\ 2 < CH 4 and E||fi 1 (pn 
p(H))||l < Ci? 4 , and therefore 



(E|| 







'q(h) 


Ph_ 




p{H) 



j!/2 < CH 2 



(65) 



where C 's are positive constants independent of e 1 but dependent on scaleless elasticity matrix 
K, scaled damping coefficient \Jec, amplitude of noise a, and slow potential V(-). 

Proof. The exact solution to the bridge dynamics is 

B n (H)q + B 12 (H)p a + j" B 12 (s)(-VV(q Q ))ds + Rq(H) 
B 21 (H)q + B 22 (H)p + jf B 22 (s)(-VV(q Q ))ds + Rp(H) 

Hence almost surely q{H) — qu = L Bi 2 (s)(—W(qo))ds. 

Since B 12 (s) — O(s) by Lemma 6.1 and E|| — VV^(go) || 1 i s bounded by Condition (3.1), one 



P(H) 



(66) 



gets 



E\\q(H) - q H \\ 2 2 < Hi E\\B 12 (s)(-VV(q ))f 2 ds 
Jo 



< H 



H 



E(||B 12 (s)|| 2 || - VV(q )\\ 2 ) 2 ds 



H [ 0(s 2 )E\\ ~VV(q )\\lds 
Jo 

0(tf 4 ) 



(67) 
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Investigation on p by applying Lemma 6.1 and Condition 3.1 gives: 



Ellfi- 1 ^)-^)^ 



E\\n-\ ( B 22 {s)ds{-VV{q {) )) + HVV(B n (H)q + B 12 (H) Po + Rq{H)))\\ 2 2 

Jo 

/ n-^B^s) - I)ds(-VV(q Q )) + fffi- 1 (VF( J B u ( J ff) go + B 12 (H) Po 
Jo 



E 



+Rq(H))-VV(q ))\\i 



< 2E[|| [ Q- 1 (B 22 (s) - I)d s {-VV(q ))\\ 2 2 + \\Hn- 1 (VV(B 11 (H)q + B 12 (H) Pa 
Jo 



+Rq(H))-VV(q ))\\ 2 2 } 



< 2[H E|jrj- 1 (B 22 (s)-/)(-VF(q ))||^s + E||i/rj- 1 (Vy(B 11 ( J ff)g 

Jo 

+B 12 (H)p + Rq(H)) - W(<?o))|||] 
Jo 

+B 12 {H)p + Rq{H)) - VV(q ))h] 

< 2H[0{H 3 ) + tfHEWn-^Bufflqo + B 12 {H) Po + Rq(H) - q Q )\\ 2 2 ] 

< 2H[0{H 3 ) + ZtfHiEWn-^BniH) - I)q \\ 2 + \\Ct~ 1 B 12 (H)p \\ 2 
+E\\n- 1 Rq(H)\\ 2 2 )) 

< 2H[0(H 3 ) + ZtfHiWn-^B^iH) - /)||lE||g ||2 + l|Bi 2 (fr)lllE|bo||l 
+E\\Rq(H)\\l)} 

< 2H[0{H 3 ) + 3L 2 H{O{H) 2 E\\q \\ 2 + 0{H) 2 E\\ Po \\ 2 + 0{H 3 ))} 
= 0{H±) 



Therefore E|| 









Ph_ 




p(H) 



\ 2 E = 0(H 4 ) independent of e 



(68) 
□ 



Lemma 6.5. Consider evolutions of different local initial conditions under the bridge dynamics: 



dqi = pi dt 

dpi = -e- 1 Kq 1 dt-VV(q 1 (0))dt-cp 1 dt + o-dW t 

dq 2 = p 2 dt 

dp 2 = -e- 1 Kq 2 dt-W(q 2 (0))dt-cp 2 dt + o-dW t 



(69) 
(70) 



Denote by L the Lipschitz coefficient of 'VV(-) (i.e. \\W(a) — W(b)\\ 2 < L\\a — b\\ 2 ), then almost 
surely 



Pl (H)-p 2 (H) 



\\e<(1 + HL)\\ 



9i(0)-9 2 (0) 

Pl(0)-P2(0) 



(71) 
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Proof. Write out the solution to the bridge dynamics in integral form: 



Pi(H) 
P2(H) 



= B(H) 
= B(H) 



9i (0)" 




Pl(0)_ 




92(0)" 









B(H-s) 



,-H 



B(H-s) 





-W(9i(0)) 


-W(g 2 (G)) 



ds 



H 



B(H-s)Y,dW s 



H 



B{H - s)Y,dW s 



(72) 



Hence almost surely 



< 



< 



< 



< 



B(H) 



p 1 (H)-p 2 (H) 

?i(0)-&(0) 
Pi(0)-pa(0) 

?i(0)-&(0) 
Pi(0)-P2(0) 

'«i(0)-&(0) 
Pi(0)-p 2 (0) 

«i(0)-&(0) 
pi(0)-pa(0) 



H 



- s) 






W(&(0))- W(&(0)) 



s (is 



W(&(0))-W(«i(0)). 


9 2 (0) ~ ft (0) 
92(0) -ft (0)' 




gi(0)-g 2 (0)' 
Pi(0)-p 2 (0) 



(73) 

□ 



Remark 6.2. If the traditional method of investigating the Lipschitz coefficient of the vector 
field is employed to evolve the separation of local initial conditions, e _1 will exhibit in the bound 
of separation. Instead we only looked at the soft part of the vector field and whence obtained a 
uniform bound. 



Theorem 3.1 (global error bound in energy norm). 

Proof. 
x(NH) 



0(H 3 ' 2 ) 
root mean square 



original dynamics 



bridge 



ejv_i(l+ifL) 
almost surely 



0(H 2 ) 

8 — — - — — xnh 

root mean square 



dynamics 



bridge 



dynamics 



1st order SIM 



Atlas of error propagation 

Let e^r = (^\\x(NH) — Xnh\\%) 1 ^ 2 - Let a and (3 be respectively the evolution of the real 
solution x((N — 1)H) and the numerical solution ^(jv_i)h by time H under the bridge dynamics 

pi. 



1G 



Then by Lemma 6.3 and 6.4 there exist constants C\ and C 2 independent of e 1 such that 



(E\\X(NH) -<5|||) 1/2 < CiiJ 3/2 

^11 1) 1/2 < C 2 H 2 (74) 



Also since ||a — /3\\e < (1 + HL)\\x((N — l)h) — £(jv-i.)/i||-E almost surely (Lemma 6.5), we 
have: 

(EH5-/3IH) 1 / 2 < (l + HL)e N -! (75) 

All in all, 

e N < (E|| a ;(7VH)-fi|||) 1 / 2 + (E|| ( i-^|||) 1 / 2 + (E||^-X J v H )|||) 1 / 2 

< {l + HL)e N _ 1 + {C 1 +C 2 )H z ' 2 

= (1 + HLfeo + ( Cl + C 2 )H^ ^^l)-! 

< (ft^^,!^^^ (76) 



Therefore letting C = ( c i+^)( e - 1 ) we have 



(E||g(T) - *fr]|i)Va < eiV < CH^ (77) 
(E||p(T) - vt\\1) 1/2 < e-^WVKhex < e-^yUhCH 1 ' 2 (78) 

□ 
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